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A numerically cheap way to obtain structural information about clusters of rare gas atoms at 
low temperatures is developed. The semiclassical frozen Gaussian imaginary time propagator is 
extended such that it can account for the mean values of all inter-atomic distances in the cluster 
and their variances. To reduce the required numerical effort an approximation for the mean values 
is developed which preserves the quality of the results offered by the semiclassical ansatz. The 
method is applied to the Are cluster. It is found that the cluster dissociates almost in one step to 
six free atoms when the temperature is increased. Precursors of the dissociation are only observable 
in the distances of the atoms via the appearance of a second isomer. The process is almost classical. 

However, the method is able to resolve small differences in the temperatures at which the dissociation 
takes place and in the mean distances of the bound configuration. 
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I. INTRODUCTION 

At very low temperatures of a few Kelvin rare gas 
atoms may assemble and form clusters due to the van 
der Waals interaction. It has been found in a variety 
of numerical studies [1-11] that these clusters exhibit a 
large number of interesting effects. In particular, struc¬ 
tural changes of the systems for an increase or decrease 
of the temperature provide insight into the behavior and 
binding mechanisms of the quantum mechanical objects 
at finite temperatures. Thus, rare gas atomic clusters 
are a topic of ongoing research. Among the rich vari¬ 
ety of their thermodynamic properties are a change from 
one packing of the atoms to another with increasing en¬ 
ergy, or phase transitions from a solid-like behavior to a 
liquid-like arrangement of the atoms [1-7]. 

For the weak van der Waals interaction between rare 
gas atoms the thermal energy corresponding to a few 
Kelvin already suffices to lead to all impacts mentioned 
above. Thus, accurate quantum mechanical computa¬ 
tions are essential to obtain reliable results. In numerical 
calculations the Boltzmann operator at inverse temper¬ 
ature P, K = exp{—PH), (due to its form also called 
imaginary time propagator) is the most important quan¬ 
tity. Its trace yields the partition function Z(/3), and 
the thermal averages of every observable O follow from 
(0)/3 = Tt{K(P)0)/Z(P). However, precise calculations 
for multi-dimensional systems are still a challenge for to¬ 
day’s numerical possibilities. For example, path-integral 
Monte Carlo methods have been used to investigate rare 
gas clusters [12-14]. At low temperatures they become 
already too expensive for a few dozen atoms, and ef¬ 
ficient but sufficiently accurate approximations are re¬ 
quired. These approximations have to reproduce quan¬ 
tum effects correctly since they will be of eminent impor¬ 
tance at low temperatures. Still many questions are open. 


* Holger.Cartarius@itpl.uni-stuttgart.de 


For example Neia and Ness [4, 6, 8, 9] might exhibit novel 
low temperature quantum effects, such as liquid-like zero 
temperature structures of Neag as compared to a solid¬ 
like structure predicted from classical mechanics [15]. 

Semiclassical methods help even with presently avail¬ 
able computation capabilities to overcome the numeri¬ 
cal drawbacks of numerically exact quantum mechani¬ 
cal algorithms and are still developed and applied in a 
wide context concerning thermodynamic properties [2- 
4, 6, 10, 11, 16-20]. Important semiclassical approxima¬ 
tions are based on the idea to restrict the quantum me¬ 
chanical wave functions to a Gaussian shape. The Gaus¬ 
sian functions are defined completely by a low number 
of parameters as the position of the center, a momen¬ 
tum of the wave packet and a width matrix. The lowest 
numerical effort is achieved with a frozen Gaussian prop¬ 
agator, i.e. a Gaussian wave function of which the width 
matrix in the exponent is fixed for all imaginary times 
(or temperatures). 

Small argon clusters are among those rare gas clusters 
which attracted large interest for a long time and again 
recently [7, 21-33]. Despite their apparent simplicity nu¬ 
merical calculations turned out to be nontrivial, and for 
the argon trimer even sophisticated path-integral Monte 
Garlo calculations could not distinguish a complete dis¬ 
sociation from structural changes within a bound sys¬ 
tem because the numerics suffered strongly from noise 
[7]. The frozen Gaussian imaginary time propagator has 
proved to provide high quality results. It could solve this 
question and identified an almost classical dissociation 
effect [16], which was confirmed [18] with a first-order 
correction to the semiclassical imaginary time propaga¬ 
tor [34, 35]. 

For the more complicated Arg cluster the specific heat 
and mean energy do not provide enough information to 
understand its whole structure. It is important to know 
which alignment of the atoms is present. A very detailed 
knowledge of the atomic structure is available via the 
inter-particle distances. To obtain their values it is very 



common to compute the radial pair correlation function 
between particles i and j at positions and Vj, respec¬ 
tively, 


Pijir) = {H\ri-rj\-r))i3 

= Tr[K{m\r^-r,\-r)]/Z{|3), ( 1 ) 


where ()/3 is the thermal average at inverse temperature 
P and K{P) = represents the imaginary time prop¬ 

agator at the same temperature [3]. This quantity pro¬ 
vides information about the distribution of the distances 
occurring at every temperature. 

It is the purpose of this paper to show that the frozen 
Gaussian method applied to atomic clusters in Ref. [16] 
can be used to determine the distances of the atoms in 
clusters in a numerically cheap and easy way. To do so, 
an extension to the semiclassical frozen Gaussian method 
will be developed. The mean values and the variances of 
the distances dij between atoms i and j can directly be 
accessed and provide a very clear information of how the 
atoms are arranged in the cluster. Since the distances 
dij of the single pairs can be calculated in parallel with 
the same Monte Garlo sampling all of these values can be 
obtained with low numerical extra cost. This provides for 
N atoms N{N — l)/2 independent values, whereas in the 
pair correlation function (1) small differences in the mean 
distances might be hidden below broad distributions. 

The example of the Arg cluster investigated in this arti¬ 
cle demonstrates the importance of the additional struc¬ 
tural information even in this relatively simple system, 
which does not undergo a structural transformation in 
the bound state. The mean energy and the specific heat 
will show that the cluster dissociates with increasing tem¬ 
perature to a system of six free atoms in one step. How¬ 
ever, precursors of the dissociation will only be observable 
in the distances between the atoms, which will indicate 
a loss of the gronnd state configuration at slightly lower 
temperatures. The variances of the distances will turn 
out to be very sensitive to the breakdown of the struc¬ 
ture. 

The further sections of this paper are organized as fol¬ 
lows. In Sec. II mean values and variances for inter¬ 
atomic distances with the frozen Gaussian method are 
introduced. For comparison the same values are intro¬ 
duced within the more flexible thawed Gaussian variant. 
Then the method is applied to the Arg cluster in Sec. 
III. After a short introduction of the system (Sec. Ill A), 
the confining sphere (Sec. IIIB), and the proper choice 
of the Gaussian width matrix (Sec. Ill G) derivatives of 
the partition function (Sec. Ill D) and the structural in¬ 
formation (Sec. HIE) are investigated. A discussion in 
Sec. IV concludes the paper. 
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II. THE FROZEN GAUSSIAN METHOD 
A. Propagator 

The method is based on a semiclassical approximation 
of the thermal operator 

K{P) = e-P^, (2) 

where /3 = I/(fcr) is the inverse temperature. The ap¬ 
proximation consists of evaluating K{P) by solving the 
Bloch equation 

= H\qo,T) (3) 

approximately for a frozen Gaussian coherent state in 
position space representation, 

{x\qo{T)) = exp|^-i[a; - q(r)]^ 

X r[x - q{T)] + {t) ■ [x - q{T)]J , ( 4 ) 

where the width matrix F is a free parameter and has to 
be adapted to the given problem as will be explained later 
for the cluster considered in this article. The propagation 
in Eq. (3) is done in imaginary time r up to the value 
T = P one is interested in. 

With Gaussian averages of the type 

/ OO 

dx^^\{x\qo{T))\'^h{x) (5) 

-OO 

the symmetrized form of the frozen Gaussian propagator 
is given by 

{x'\Kfg{t)\x) = det(r)exp ^-■^Tr(r)T^ 

X ^ det ^2 [1 — exp(— ^.^rr)]”^^ 

X exp i[a;' — a;]"'"r[tanh(fi.^rT/2)]“^[a;' — a;]^ 

-[x-q{T/2)]^T[x-q{T/2)]^, (6) 

where x = (a;' -|- x)/2. The partition function is simply 
given by the trace [35] 

Zfg{t) = Tr[A:FG(T)] = \/det(r)exp ^-■^Tr(r)Tj 
X J det ^[1 — exp(— 
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/■°° / /•^/2 \ 

The whole dynamical information is contained in the 
imaginary time propagation of the variable q{T) and is 
governed by the SN coupled equations of motion 

T = -r-i(Vf^(q(T))), (8) 


directions are equivalent the best choice are multiples of 
the unity matrix Isxs, such that two scalar coefficients 
Di and D 2 have to be chosen, i.e. 

Di = Di1sx3, D 2 = I?2l3x3- (11) 

The form (6) of the propagator is represented in Carte¬ 
sian coordinates and is the most efficient for the numer¬ 
ical evaluation. Consequently the Cartesian representa¬ 
tion 


which are relatively simple and can be integrated with 
a standard integrator for ordinary differential equations. 
The remaining numerical task is a single position space 
integration for the initial positions q(0) in Eqs. (6) or (7), 
which is done with a Monte-Carlo integration. 

As shown previously [16, 18] a reasonable choice of the 
width matrix is crucial for the quality of the semiclassical 
method. However, highly precise values can be obtained 
already with a very simple structure. Since all particles 
are identical and thus the pairwise interactions are the 
same for all combinations only the center of mass motion 
has to be distinguished. In center of mass coordinates. 



i=l 


i = l...,Af-l, (9b) 


good semiclassical estimates are obtained with the matrix 


( D AD AD 
AD D AD 

V 

D = {Di + {N -1)D2)/N, (12b) 

AD = {D^ - D 2 )/N (12c) 

of the matrix (10) is used. 

B. Thermal averages of structural information 

Simple thermal averages requiring only the partition 
function Z(/3) are the mean energy E = YT'^dln Z /dT 
and the specific heat C = dE/dT. However, K{I3) pro¬ 
vides access to the thermal average of any observable O 
via 



-(PG) ^ Tr(KFG(/3)0) 

Zfg{P) 


(13) 


r 


cmc 


(Di 0 
0 D 2 

V : 


( 10 ) 


and the 3x3 diagonal matrices Di and £>2 describing 
the three spacial directions of the center of mass and the 
relative coordinates, respectively. Since also the spacial 

_I 


which is exploited in this article to gain access to the 
structural information. A well suited property is the dis¬ 
tance between two atoms, i.e. 

O = dij = \Xi I ■ (^^) 

With the frozen Gaussian propagator (6) this leads to 
the expression 


JIFG) ^ _1_ 

ZfgW) 


Tr{KFG{P) \xi - Xj\) = 


1 


Zfg{P) 

13/2 


det(r)exp ■^Tr(r)/3^ ^det ^2 [1 — exp(—fi^F/?)] 


/-oo (27r)3^ 


exp(-2^ dT(H(g(T)))^ / e^p(^-[x - q{P/2)fT[x - q{P/2)]^\x, - Xj\, (15) 


in which an explicit integration over the 3N position vari¬ 
ables X remains in addition to the evaluation of the parti¬ 
tion function (7). As mentioned previously [18] it is very 
important to reduce the numerical effort as much as pos¬ 
sible for many-particle systems. In particular, the posi¬ 
tion space integrations require an expensive Monte Carlo 
sampling in a high-dimensional configuration space. 


For usual applications a numerical evaluation of the x 
integration can be avoided in a reasonable approxima¬ 
tion. This can be seen with the variable y = x — q(/3/2), 
which transforms the x integral in (15) to 

/ = y'“ dy3^exp(^-yTry^ 


X Ivi - Vj + g*(/3/2) - qj{P/2)\. (16) 
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The widths of the atom’s wave functions contribute only 
at low temperatures significantly to the distance. As will 
be seen, in practical applications a very narrow Gaussian 
centers all values yt strongly around zero, i.e. Xi is al¬ 
most identical with qi{P/2) for a nonvanishing Gaussian 
weight. Thus, the integral (16) is calculated for the case 
IVi — Vjl ^ \Qi — Qj\- With the expansion 


\y^ 


Vj -b q* - QjI « 

+ \{y^-y3f 


1 [{yi -yj)- - Qj)f 

9 


the integral evaluates to 


r3N 


1 = 


det(r) 


\q,{l3/2)-q,{/3/2)\ 


Tr(r)-i 


. r-.i - r".i 

JJ ^3 


- 

} 


Q\qi{P/2) - q0/2)\ 


(18) 


where is the 3x3 submatrix of F at the rows and 
columns representing particles i and j. In total 


% 


1 


( 3 ) = \x^ - Xj\) « \/det(r)exp ^-^Tr(r)/3^ ^det ([1 - exp(-?i2r/3)] 


■^fg(/3) 

/.oo 


(27r)3^/2 


r/5/2 


exp I —2 


dT{V{q{T))) 


\q3P/2)-q,iP/2)\ 


6|q*(/3/2) - qj(/3/2)| 


(19) 


is obtained. 

The first term in equation (19), oc \qi — qj\, reflects 
the core of the semiclassical approximation, in which the 
positions of the atoms are given by the centers qi of the 
Gaussian wave packets (4). It corresponds to 

0 = \q^-qJ\■ (20) 

The second term contains a correction due to the finite 
width of an atom’s wave packet. It is completely suffi¬ 
cient to include this lowest-order term, of which the x 
integration could be done analytically with a simple re¬ 
sult, thus reducing the numerical effort drastically. For 
the frozen Gaussian method any higher terms beyond 


those included in the approximation (19) for the mean 
distances are of lower interest. From the physical point 
of view it is expected that the width of the atom’s wave 
function only plays a role at very low temperatures at 
which the structural configuration is unambiguously in a 
highly symmetric ground state configuration. Indeed, as 
will be seen in the results already the correction term in 
the approximation (19) is very small. 

To estimate the quality and validity of these mean val¬ 
ues additionally the variances of the distance distribu¬ 
tions are calculated. With the operator 

O = Vij = {xi - Xjf - (21) 

and the integral 


^ Tr(A:FG(/3) [xt - Xj]'^) = \/det(r)exp ^-■^Tr(r)/3^ ^det (^[1 - exp{-h?ri3)] 


^FG iP) 


dq 


3N 


(27r)3^/2 


ZfgW) 

r/3/2 


expl —2 


dT(F(q(T))) 


(q.(/3/2) - q,(/3/2))2 + + F"/ - F".i - F"/) 


( 22 ) 


the standard deviations aij 
are obtained. 

‘'3 



of the distances 


C. Sorted distances 


For the Are cluster there are 15 possible combinations 
i and j, and thus 15 distances. The clusters are oriented 
arbitrarily in the simulation. The numbers i and j of the 
atoms have no meaning for the true configuration, and 
thus are not appropriate quantities to define the pairwise 
distances. The average of all calculations simply results 
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in identical values for all dij, which correspond to the 
mean value of all 15 atom-atom distances in a certain 
configuration. To obtain a meaningful quantity the dis¬ 
tances are sorted according to their size, 

di < d2 < ■ ■ ■ < di5, (23) 

and the thermal average of these size-ordered distances is 
determined, i.e. the thermal mean values of the smallest 
distance, the second smallest, and so forth are obtained. 
These values can be compared with the expectations of 
geometrical configurations. In an experiment the single 
distances are accessible [36] and can in a given sample 
be sorted the same way. Alternatively, results from this 
calculation can be used to determine the distance of the 
atoms with a well-grounded assumption about the con¬ 
figuration [37]. 


D. Comparison with the thawed Gaussian 
propagator 

The frozen Gaussian method has proved to provide 
good results for thermodynamic quantities. We want to 
know whether or not this is also true for the widths cal¬ 
culated in this article. Thus, the structural information 
of the frozen Gaussian method is compared with that of 
a more flexible thawed Gaussian ansatz. It is based on 
a time-dependent width matrix G(t), which adapts it¬ 
self to the given temperature. With the restriction to 
Gaussian wave packets the thawed Gaussian variant is 
usually the most accurate approximation. The variable 
width matrix adds an additional freedom in the param¬ 
eters. This is reflected in the quality of the results as 
has clearly been demonstrated for a double well potential 
[38]. For a large number of degrees of freedom it suffers, 
however, from the higher numerical costs. The single¬ 
particle ansatz of Frantsuzov et al. [3] avoids these diffi¬ 
culties by reducing the matrix G(r) to a block-diagonal 
structure, where 3x3 matrices representing the three spa- 
cial coordinates of one particle are the only non-vanishing 
matrix elements. In the case of six atoms this reduction 
is not required and there is no need to ignore the inter¬ 
particle correlations. 

The thawed Gaussian propagator used for compari¬ 
son with the frozen Gaussian method is the time evolved 
Gaussian approximation (TEGA) suggested by Frantsu¬ 
zov et al. [3, 39] with a full width matrix G. It is based 
on the solution of the Bloch equation (3) with the coher¬ 
ent state 

{x\g{q{T),GiT))) = (7r3^1detG(T)l) 

X exp - q{r)VG{ t)~^[x - ^(t)]^ . (24) 

The resulting symmetrized propagator reads 


X exp - q(r/2)]^G(T/2) ^[a? - q(T/2)]^ 

X exp (^-^[®' - q{r/‘^)VG{ t/ 2)-^[x' - q(T/ 2 )]^ 

(25) 

with the time-dependent width matrix G{t). In imag¬ 
inary time T the equations of motion for the Gaussian 
parameters G, q, and 7 are 

|-G(r) = -G(r)(VVTE(q(r)))G(T) + hH, (26a) 

^q(r) =-G(r)(VE(q(r))), (26b) 

^ 7 (r) = -^Tr [(VV'^I/(q(T)))G(r)] - (E(q(T))), 

(26c) 

which have to be integrated from r = 0 to larger times 
with the initial conditions 


q(T « 0) = qo, (27a) 

G{TK.{)) = h^lT, (27b) 

7 (t « 0) = -E(qo)T. (27c) 


In all expressions (...) represents Gaussian averaged 
quantities of the form (5) with the wave packet (24), and 
1 is the iN X 3A^-dimensional identity matrix. The rel¬ 
evant quantities are the partition function 


^ f dq^^ exp[27(T/2)] 

™ J (2v/^)3^ Vdet[G(r/2)] 


(28) 


and the mean value of the distances in the same approx¬ 
imation as for the frozen Gaussian method. 




dq' 


3N 


exp[ 27 (/ 3 / 2 )] 


Ztg(/3) J (20F)3^ v'det[G(/3/2)] 
lq,(/3/2)-q,(/3/2)l 


-f 


Tr(G,,(/3/2) 


.G,7(/3/2)-G.,(/3/2)-G,.(/3/2)) 


6lq.(/3/2)-q,(/3/2)l 


(29) 


The thawed Gaussian approximation allows for an 
additional important information. Its temperature- 
dependent width matrix G(t) provides easier access to 
the width of the wave function, which influences the vari¬ 
ances of the distances. The quantum mechanical part 
of the variances, i.e. that originating from the spread of 
the wave function, is expected to increase at lower tem¬ 
peratures. For a frozen Gaussian this can be described 
correctly if the constant matrix F is optimized for ev¬ 
ery single temperature. In the thawed Gaussian case the 
variances read 


{x\Ktg iT)\x') 


dq3^ exp[27(T/2)] 
(2^)3^ det[G(T/2)] 


1 f dq3^ exp[ 27 (/ 3 / 2 )] 

ZtgW) J (20F)3^ ^det[G(/3/2)] 
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(q.(/3/2) - q0l2)f + -Tr(G.i(/3/2) + G,,(/3/2) 


- G.,(/3/2) - G,,(/3/2)) 




and follow directly from the imaginary time evolution of 
q{T) and G{t). We are mainly interested in the quantum 
mechanical part of the variances, viz. 


-(TG.qm).^, _ 1 f dq^^ exp[27(/3/2)] 

^ Ztg{P) J (2v^) 3^ v/'det[G(/3/2)] 

X iTr(G„Ci/2) + aijiP/2) - atiWft) - G„(^i/2)). 

(31) 


III. STRUCTURAL INFORMATION ABOUT 
THE Are CLUSTER 


A. Representation of the system 

The argon cluster consists of 6 atoms, where the Hamil¬ 
tonian in mass scaled coordinates reads 

iL = -y^A, + ^U(ry) (32) 

i—1 j<i 

with the Laplacian of particle i. The two-body po¬ 
tential V{rij) of Argon is still a very challenging task. 
One of the best analytic expressions at hand is a fit to 
experimental results by Aziz and Slaman [40] of which 
an adoption to a Morse potential [25] is used, 

V{rij) = D (exp [-2a{rij - i?e)] - 2 exp [-a{rij - Re)]) 

(33) 

with the parameters D = 99.00cm a = 1.717A, and 
Re = 3.757 A in consistence with previous studies of the 
Argon trimer [7, 16, 18]. 

The numerical efficiency of the frozen Gaussian 
method is increased with an expansion of the potential 
in terms of Gaussians, viz. 

Vi\ri-rj\) = , r,j = In - rj], (34) 

p 

This procedure was suggested by Frantsuzov et al. [3] 
and has successfully been applied [3, 16, 18]. In the form 
(34) Gaussian integrals of the potential or its derivatives 
can be done analytically. The required parameters for 
a fit to three Gaussians are listed in Table I and were 
previously obtained in Ref. [16]. 


B. Confining sphere 

An additional potential is usually introduced to con¬ 
verge the numerical q integration. All particles are con¬ 
fined within a sphere around the center of mass Rem by 


TABLE I. Parameters of the Argon-Argon interaction poten¬ 
tial expressed in terms of Gaussians according to Eq. (34) 
[16]. 


p 

Cp [cm“^] 

[A""] 

1 

3.296 X 10'" 

0.6551 

2 

-1.279 X 10® 

0.1616 

3 

-9.946 X 10® 

6.0600 


the condition \q — ilcm| < Rc, where Rc is the confining 
radius. This can be achieved with the steep potential 

2=1 ^ 

added to the Hamiltonian (32) or, as in our study, by a 
restriction of the volume for the q integration. 

Of course, an additional potential influences the re¬ 
sults and can crucially change the behavior of the cluster 
[2, 21]. If only bound configurations are investigated Rc 
is usually chosen such that the bound configurations are 
not affected, i.e. Rc is larger than the extension of the 
bound cluster. However, we are interested also in the 
dissociation process for which the choice of Rc is non¬ 
trivial [16, 18]. A larger radius Rc always allows for a 
dissociation at lower temperatures. In principle it has to 
be adopted to the physical conditions as, e.g. the pres¬ 
sure. We are interested in the qualitative behavior at the 
dissociation and it was checked carefully that the choice 
i?c = 35 A does not influence the qualitative change of 
the relevant observables, i.e. the mean energy, the spe¬ 
cific heat, the mean values of the inter-atomic distances 
and their variances. In particular, it was assured that 
the case of a completely dissociated cluster is present for 
temperatures above 40 K and the form of the dissocia¬ 
tion process is not altered. The value of the confining 
radius Rc, i.e. the pressure in physical terms, affects the 
temperature at which the dissociation occurs. 


C. Choice of the width matrix 

While in a thawed Gaussian calculation the initial con¬ 
dition for the width matrix (27b) is defined, the constant 
matrix F of its frozen Gaussian counterpart has to be 
chosen carefully. It is a free parameter of the system. 
It is not trivial to find a good choice of F. However, as 
was mentioned above, the structure (12) with the 3x3 
submatrices (11) is well suited. Thus, only the two pa¬ 
rameters Di and D 2 need to be chosen. 

A detailed investigation of the Ara cluster revealed 
that there is a reliable and simple method to find the 
best choice for the inter-particle width parameter D 2 [16]. 
Propagating the partition function (7) to long imaginary 
times /3 —>■ 00 , i.e. T —>■ 0, one can extract the thermo¬ 
dynamic mean energy to correspond to the ground state 
energy Eq . The parameter D 2 providing the lowest value 
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for i?o has shown to lead to the best agreement with nu¬ 
merically exact calculations and the more flexible thawed 
Gaussian approximation. This result is almost indepen¬ 
dent of the temperature at which the partition function, 
mean energy or specific heat of all methods are compared. 
Furthermore, calculating the first-order correction to the 
frozen Gaussian propagator (6) showed that this choice 
also requires the smallest correction. Thus, the simple 
minimization of the ground state energy gives us a reli¬ 
able way of determining 02- For the Are cluster in this 
article it was found that £>2 = 32 A is the best choice. 

The center of mass is free and it can exactly be de¬ 
scribed by a Gaussian in the limit Z3i —^ 0. This means 
the value should be as small as possible. For the numer¬ 
ical evaluation one needs a finite value. It is known that 
£>1 = 0.1 A"' is small enough [16]. The results cannot 
be distinguished from those of even lower values for Di. 

D. Mean energy and specific heat 

To compare the Arg cluster with the trimer consid¬ 
ered in [16, 18] the mean energy and the specific heat 
are studied first. They are shown in Fig. 1 for the two- 
parameter (£>1 and D 2 ) frozen Gaussian propagator and 
the derivatives of the classical partition function 

The observations are very similar to those obtained for 
trimer. At very low temperatures the classical calcula¬ 
tion tends to the potential minimum. At T = 1K a 
mean energy of E = — 1216cm“^ is found. The frozen 
Gaussian results are at this temperature already in a 
very flat regime, in which the mean energy is almost 
independent of the temperature and approximates the 
quantum mechanical ground state energy. The method 
leads to Eq fv —1015cm“^. In the inset of Fig. 1 (a) a 
comparison with the more flexible fully-coupled thawed 
Gaussian propagator is shown. It leads to a value of 
Eq fv —1040 cm“^, i.e. the difference of the ground state’s 
binding energy is only 2.4%. Thus, one may conclude 
that also for the larger Arg cluster the quality of the 
frozen Gaussian propagator is acceptable in comparison 
with the numerically more expensive thawed Gaussian 
variant even for the low-temperature limit. 

The dissociation appears in the mean energy as a step. 
The energy raises almost directly to that of six free parti¬ 
cles. This indicates a dissociation of all atoms at once as 
was observed for the trimer. The same information can 
be gained from the specific heat, which is shown around 
the dissociation in Fig. 1 (b). One broad peak confirms 
that the dissociation occurs in one step. The classical 
calculation shows a transition at a slightly lower temper¬ 
ature, and the difference between the two maxima in the 
specific heat is approximately 0.5 K, which is lower than 
for the trimer, where a difference of 1.5 K was observed. 



r[K] 



r[K] 


FIG. 1. (Color online) (a) Mean energies of the Arg clus¬ 
ter calculated with the two-parameter frozen Gaussian prop¬ 
agator (2P-FG) and its classical counterpart. For the high- 
temperature limit both results agree well. At low tempera¬ 
tures the classical mean energy tends to the potential mini¬ 
mum and the frozen Gaussian propagator approximates the 
quantum mechanical ground state. The inset shows a com¬ 
parison of the 2P-FG method with the fully coupled thawed 
Gaussian method (FC-TG). The energies differ by a few per¬ 
cent. (b) Specific heat around the dissociation, which seems 
to happen in one step. 


Gertainly the difference can again be related to the zero 
point energy in quantum mechanics. This is larger for 
six atoms than for three and one could expect that also 
the temperature difference is larger. However, one has 
to keep in mind that this energy has to be distributed 
among a larger number of atoms during the dissociation. 
Since aside from the small shift in the temperature the 
dissociation process is almost identical in the classical 
and the frozen Gaussian calculation one may conclude 
that it is a purely classical phenomenon. 


E. Structural information for low temperatures 
and for the dissociation 

In addition to the information of the simple derivatives 
of the partition function the mean distances are studied. 
Since so far the dissociation seems to be purely classical 
it is interesting to also compare the structural informa¬ 
tion with the classical one. To do so, the classical mean 
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FIG. 2. (Color online) (a) Comparison of the mean values of 
all 15 distances calculated with the frozen Gaussian method 
(2P-FG) and a classical calculation at temperatures T < 20. 
The distances appear in groups. Three distances converge 
for T —>■ 0 to values above 5 A, whereas the remaining 12 
are below 4 A. (b) A comparison of the frozen Gaussian and 
thawed Gaussian method (FC-TG) shows that the mean dis¬ 
tances agree very well. 


distances 


^classical) 

% 


1 / kr r 
Zci \27rft2y J 


.-/3T(q) 


- qj\dq 


.3N 


are added in the calculations below. 


(37) 


1. Structure at low temperatures 

The mean distances obtained for temperatures below 
20 K, i.e. significantly below the dissociation process, are 
shown in Fig. 2, where first the frozen Gaussian method 
(2P-FG) is compared with the classical calculation in 
Fig. 2 (a), and then the fully-coupled thawed Gaussian 
approximation is added in Fig. 2 (b). The most strik¬ 
ing observation is that the distances appear in the low- 
temperature limit in two groups. A group of three “long” 
distances, of which the values are always above 5 A for 
T —>■ 0, and a second group of the 12 remaining “short” 
distances, which converge to a value below 4 A, exist. 
This already gives a clear answer to the question about 
the ground state configuration of Arg. It is consistent 
with the distances in an octahedron, or in other words, 
the atoms are located at the centers of the surfaces of its 
dual polyhedron, viz. the cube. Twelve short distances ds 


from the atoms on neighboring surfaces and three longer 
distances di between the atoms on opposite surfaces are 
expected. The ratio of the distances is supposed to be 
di = y/2ds, which is fulfilled excellently in both the clas¬ 
sical and semiclassical calculations. 

In the classical case the atoms seek directly the po¬ 
tential minima, whereas in the quantum case always a 
wave function with a finite width leading automatically 
to larger mean distances is present. Furthermore, in the 
classical case the fixed octahedron configuration is only 
observable for T —>• 0. This is a consequence of the fact 
that classically every nonvanishing energy allows for a 
thermal excitation. In contrast to this there should be no 
excitation possible if kT is clearly below the energy differ¬ 
ence between the ground state and the first excited state 
in the quantum mechanical case. This is also reflected 
in the mean distances. For temperatures below T Ri 3 K 
no differences between the distances in one group are ob¬ 
served, and the distances do not change for even lower 
temperatures. This indicates that the cluster is already 
in the ground state configuration. 

A comparison of the frozen Gaussian method with the 
fully-coupled thawed Gaussian propagator reveals that 
the distances agree very well. This is in particular true for 
all larger distances. Also the low temperature limit shows 
an excellent agreement. The edge length of the cube con¬ 
taining the octahedron is di = 5.39 A in the frozen Gaus¬ 
sian calculation and di = 5.35 A in the thawed Gaussian 
approximation. The difference is below 1 %, and thus 
even smaller than that of the mean energy. Obviously 
the structural information of the frozen Gaussian method 
is less affected by the constant Gaussian width approxi¬ 
mation. 

With the data of Fig. 2 we are also able to estimate 
the quality of the approximation (17), in which the power 
series expansion of the distances was introduced. The 
first-order term retained in the expansion is of the size 

TrfT'^ -h TT-i - r~-i - 

G\qi-qj\ 

where for the width matrix (12) Tr(r))^ -|- — 

r“-^)/6 = D 2 ^ = (32A = 0.031 A^ is obtained. 

With the knowledge that the typical distances \qi—qj \ are 
even in the bound phase of the order of a few Angstroms 
this correction can be estimated to be always less than 
1 % of the leading order. Hence, it has at most the same 
size as the difference between the two semiclassical prop¬ 
agators. Higher orders in the series expansion (17) would 
lead to even smaller corrections, which do not need to be 
taken into account since they are below the error of the 
semiclassical approximation. 


) 


2. Dissociation to six free atoms 

The distances around the dissociation are shown in Fig. 
3, where the classical and 2P-FG results are compared. 
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FIG. 3. (Color online) (a) Mean values of all 15 distances 
calculated with the frozen Gaussian method (2P-FG) around 
the dissociation. Two of the smaller distances join the group 
of the three larger distances starting at T « 23 K, but then 
the dissociation occurs at once, (b) The same behavior is 
observed for the classical calculation. 


Since the frozen Gaussian approximation is known to pro¬ 
vide good results at these temperatures (cf. Ref. [18]), a 
comparison with the thawed Gaussian propagator does 
not give any new information. Figure 3 confirms the 
finding of the consideration of the mean energy and the 
specific heat in Fig. 1. The dissociation effect is clas¬ 
sical. The semiclassical approximation of the quantnm 
mechanical propagator and the purely classical calcula¬ 
tion lead to the same behavior. Apart from a small shift 
in temperature both diagrams agree very well. 

The calculation of the distances gives additional insight 
into the dissociation process. For temperatures T ^ 23 K 
the groups of short and long distances are unchanged. 
Above this temperature two of the twelve distances ds 
are separating from the others and join the three longer 
distances di. Two new groups with five and ten distances 
start to from. For Arg a second isomer in the form of a 
tri-tetrahedron is known to contribute at increasing tem¬ 
peratures [33]. This would exactly agree with a grouping 
of five longer and ten shorter distances and is also visible 
in Fig. 4, in which the distances of all six atoms from 
the center of mass are shown for a classical calculation. 
At low temperatures all distances have almost the same 
size, which agrees with a pure octahedron configuration. 
With increasing temperature one distance grows and in¬ 
dicates a coexistence of octahedron and tri-tetrahedron 
configurations. Signatures of further arrangements of the 



FIG. 4. (Color online) Distances of all six atoms from the 
center of mass in a classical calculation for Rc = 35 A. Shown 
is a thermal average of size-ordered distances as explained in 
Sec. IIG. With increasing temperature one distance becomes 
larger than the others before all raise drastically. 


atoms are not found. 

However, this rearrangement of the atoms remains in 
its beginnings. A new structure cannot completely arise 
since the whole process does not finish before all distances 
raise drastically and indicate with this increase the dis¬ 
sociation of the cluster. Comparisons with calculations, 
in which due to a smaller value of Rc no dissociation is 
allowed, show that this effect only appears in connec¬ 
tion with the dissociation. Thus, the separation of the 
two distances is more a precursor of the total destruction 
of the cluster. The dissociation happens then at once. 
After the dissociation the distances obtain new almost 
constant values which correspond to the distribution of 
atoms moving freely within the confining sphere. 


3. Variances of the distances 

To learn more about the actual distribution of the dis¬ 
tances their standard deviations are shown in Figs. 5 (a) 
and (b) for the classical and the 2P-FG method. The 
most significant feature is the drastic increase of the stan¬ 
dard deviations around the transition. At this temper¬ 
ature range parts of the simulated clusters are still in a 
bound configuration whereas others are already dissoci¬ 
ated. The effect is similar and of the same size for the 
classical and the semiclassical calculation. For tempera¬ 
tures above the dissociation the standard deviations are 
almost the same for all distances, which is also expected 
for six free atoms. 

Of more interest is the behavior of the standard devi¬ 
ations below the dissociation. The longer distances are 
expected to show more fluctuations. Additionally two of 
the shorter distances join the group of the three longer 
distances for increasing temperatures as a consequence of 
contributions from two isomers. It can also be expected 
that these two show higher standard deviations than the 
short distances since the separation of the two distances 
does not happen abruptly at one temperature as can be 
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FIG. 5. (Color online) (a) Standard deviations aij of the dis¬ 
tances in the frozen Gaussian approximation (2P-FG) and (b) 
of the classical calculation. They increase drastically around 
the dissociation, (c) Quantum mechanical part (31) of the 
variances for the fully coupled thawed Gaussian propagator. 
For low temperatures it increases due to an increasing width 
of the Gaussians. 


seen in Fig. 3. Consequently, below the dissociation the 
standard deviations are found to form two groups. One 
group combines the standard deviations of five distances, 
i.e. the three longer ones and the two joining them. The 
other group consists of the standard deviations of the 10 
short distances which stay together up to the tempera¬ 
ture of the dissociation. 

In Figs. 5 (a) and (b) it seems that the standard de¬ 
viations always decrease for lower temperatures. This 
is definitely expected for the classical calculation. The 
frozen Gaussian approximation can, since the width of 
the wave function is determined by the constant values of 
r, not reflect the quantum mechanical expectation that 
the wave function stretches at lower temperatures. To 
cover also this effect the quantum mechanical part ac¬ 
cording to equation (31) is plotted in Fig. 5 (c). Only 
the variance of one of the distances is shown since this 
part is almost identical for all of them in the highly sym¬ 
metric situation of a monoatomic cluster. At T 20K 
the variance has approximately the same size as that 
following from the frozen Gaussian method = 


Tr(r))i -b r-/ - r-.i - r-/)/2 = = 0.093 Al in 

particular, these quantum mechanical parts of the vari¬ 
ances are considerably lower than other contributions in 
equations (22) and (30). Figure 5 (a) would not change 
with the thawed Gaussian propagator. Only for tem¬ 
peratures T ^ 10 K the extension of the wave function 
becomes important for a measurement of the distances. 
For temperatures in the range of the dissociation the in¬ 
formation from the frozen Gaussian method is completely 
sufficient. 

All calculations shown in this section could be imple¬ 
mented and performed easily on a NVIDIA Tesla G2070 
GPU. On this architecture a converged result for low 
temperatures is obtained in less than 20 hours. The 
most critical part is that around the dissociation, where 
a very detailed sampling for bound configurations has to 
be done alongside an inclusion of large distances allowing 
for an unbound cluster. This can require an increase of 
the sampling points by a factor of 10. 


IV. SUMMARY AND OUTLOOK 


In this article it was shown that structural informa¬ 
tion about a cluster of atoms can be obtained with the 
frozen Gaussian semiclassical method in a numerically 
cheap way. The evaluation of the corresponding integrals 
can be done in parallel to that of the partition function. 
With this method the full information about all the dis¬ 
tances of all combinations of the atoms can be obtained. 
A comparison with the more flexible thawed Gaussian 
propagator revealed that the quality of the distances is 
on the same level as that of the mean energy or the spe¬ 
cific heat, or even better. To avoid inefficient numerical 
computations of a position space integral an approxima¬ 
tion for the distances was introduced. It was possible, 
however, to show that this approximation does not re¬ 
duce the quality of the results below that obtained in the 
semiclassical approximation of the propagator. 

On the physical side it was found that with increasing 
temperatures the Arg cluster undergoes an almost direct 
transition to six free atoms. However, it shows precur¬ 
sors in the distances. At temperatures slightly below the 
dissociation a reordering of the atoms starts, in which 
contributions from a second isomer, viz. a tri-tetrahedron 
[33], appear, but then vanish in the increasing distances 
at the dissociation. The dissociation is a purely classi¬ 
cal effect. The semiclassical approximation shows exactly 
the same behavior with just a small shift in the tempera¬ 
ture of 0.5 K. Around the dissociation the standard devi¬ 
ations of the distances are almost completely determined 
by classical contributions. Only for lower temperatures 
the extension of the wave functions becomes important 
as was seen in a thawed Gaussian approximation. 

At low temperatures the cluster assumes the shape of 
an octahedron, where the longer distance between the 
atoms is di = 5.4 A, and the shorter has the value dg = 
dily/2 = 3.8 A. Glassically the fixed configuration is only 
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obtained in the limit T —>■ 0 whereas in the quantum 
mechanical case the ground state configuration is present 
for all temperatures T ^ 3K. 

The frozen Gaussian method has proved to provide re¬ 
liable results for quantum mechanical calculations. There 
is a large number of investigations which can be done 
with it. In particnlar, the results for Ara in [16] and for 
Are in this work indicate that the confinement to very 
small spheres usually applied in the calculation of the 
partition function and values deduced from it [2-4, 41] 
only sample bound cluster configurations. This is phys¬ 
ically realized at high pressures. If one is interested in 
lower pressures, at which a dissociation is allowed, this is 
too restrictive to fnlly understand the low-temperature 
behavior of the clusters. The dissociation can set in be¬ 
fore structural changes or a melting can be observed. To 
take this into account it is necessary to advance the in¬ 
vestigations done here to clusters with higher numbers of 
atoms. In particular, the cases of Aris [42-44], Neis [3] 
or Ness [4] are of special interest since they showed inter¬ 
esting structural transformations in the non-dissociated 
cases. Whereas the stronger quantum effects in the com¬ 
pletely bound case are well covered by a variable width 
matrix the numerically cheaper frozen Gaussian method 
has advantages in the numerically more challenging case 


of the dissociation requiring a sampling of bound and 
unbound configurations of the atoms. The exactness of 
both methods can then be monitored and (if necessary) 
improved with the the series expansion of the imaginary 
time propagator [18, 34, 35]. 

The importance of the series expansion is not restricted 
to the dissociation. Most effects in rare gas clusters such 
as structural transformations or dissociations appear at 
such low temperatures that it is necessary to analyze 
whether the semiclassical approximations used in the cal¬ 
culations correctly reproduce the true quantum mechan¬ 
ical behavior. An important example will be Nesg, for 
which strong differences are found between the approxi¬ 
mate quantum computations and a purely classical the¬ 
ory [6]. 
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